where<-Sys.info()["sysname"]
if (where=="Darwin"){
  setwd("~/evo-dispersal/People/Louis/HZDisp1D/")
  #system(paste("R CMD SHLIB Density1D.c"))
  system(paste("R CMD SHLIB PointMetrics1D.c"))
}



#system(paste("R CMD SHLIB PointMetrics1D.c"))

dyn.load("PointMetrics1D.so")

#is.loaded("sum_metrics")

### iterates pop over ngens and collects resulting dvec #
### then re-initialises population, using dvec, and runs with and without dispersal evolution #

mother<-function(n, spX, K, lambda, B, strength, ngens, Dmean, h2H, h2D, VPH, VPD){
  pop<-init.inds(n=n, spX=spX, Dmean=Dmean,
                 h2H=h2H, h2D=h2D, VPH=VPH, VPD=VPD)
  ztrace<-c()
  for(ii in 1:ngens)
    {
    print(paste(ii, nrow(pop), sep=" "))
    pop<-repro.disp(popmatrix=pop, K=K, lambda=lambda, B=B, strength=strength, 
                    h2H=h2H, h2D=h2D, VPH=VPH, VPD=VPD, spX=spX)
    if (is.null(pop)) break
    }
  out<-list(parameters=list(ngens=ii, K=K, lambda=lambda, B=B, strength=strength, n=n, spX=spX,
                            ngens=ngens, Dmean=Dmean, h2H=h2H, h2D=h2D, VPH=VPH,
                            VPD=VPD), pop=pop, ztrace=ztrace)
  out
}
##### #
# returns the bounded scale (0-1) from the continuous scale
inv.logit<-function(cvec){
	exp(cvec)/(1+exp(cvec))
}

#returns the continuous scale from the bounded 
#unnecessary, here only for a sense of symmetry
logit<-function(pvec){
	log(pvec/(1-pvec))
}

### Initialises a matrix with n individuals randomly placed at a location X on a -spX to spX linear space #

init.inds<-function(n, spX, Dmean, h2H, h2D, VPH, VPD)
  {
  X<-runif(n, -spX, spX) #position
  D<-rnorm(n, Dmean, sqrt(h2D*VPD)) # individual sd
  HI<-HI(X)
  DP<-D+rnorm(n, 0, sqrt((1-h2D)*VPD))  
  cbind(X, D, HI, DP)		
  }
##### #


### init HI #

HI<-function(X)
{
  Hi<-vector("numeric",length(X))
  Hi[X<0]<--5
  Hi[X>0]<-5
  Hi
}
##### #


### Fitness based on hybrid index (remplace habitat variance) #

wH<-function(HI, strength, B)
  {
  1-strength*(4*HI*(1-HI))^B
  }
##### #

wAbs<-function(popmatrix, dens, strength, B, lambda, K){
	hyb.w<-wH(inv.logit(popmatrix[,"HI"]), strength, B)
	hyb.w*hass.comm(lambda, K, dens)
}

### #

repro.disp<-function(popmatrix, K, lambda, B, strength, h2H, h2D, VPHI, VPD, spX){
  SDVED<-sqrt(VE(h2D, VPD)) # define environmental variances
  SDVEHI<-0
  # calculate traits through space
  mets<-ca.metrics(popmatrix, spX)
  
  # reproduction
  Eoffs<-wAbs(popmatrix, mets[,"Nx"], strength, B, lambda, K) #expected offspring
  offs<-rpois(nrow(popmatrix), Eoffs) # realised offspring
  if (sum(offs)==0) return(NULL) #catch extinctions
  inds<-1:length(offs) #vector indexes for offspring
  inds<-rep(inds, times=offs)
  popmatrix<-popmatrix[inds,] #replace parents with offspring. Offspring inherit location and P from parents
  mets<-mets[inds, ] #expand mets matrix as well
  #mets<-matrix(mets[inds,], ncol=ncol(mets), dimnames=list(NULL, colnames(mets))) wtf?
  if (class(popmatrix)=="numeric") popmatrix<-matrix(popmatrix, ncol=4, dimnames=list(NULL, c("X", "D", "HI", "DP"))) #cases where only one individual left
  
  # calculate new G and P
    SDD<-mets[,"SDD"]
    SDH<-mets[,"SDH"]

	popmatrix[,"D"]<-rnorm(nrow(popmatrix), mets[,"MeanD"], SDD)
	popmatrix[,"DP"]<-popmatrix[,"D"]+rnorm(nrow(popmatrix), 0, SDVED)
	popmatrix[,"HI"]<-rnorm(nrow(popmatrix), mets[,"MeanH"], SDH)
  
  # Disperse offspring
  disp<-rnorm(nrow(popmatrix), mean=popmatrix[,"X"], sd=exp(popmatrix[,"D"])) #disperse offspring
  disp<-ifelse(disp<(-spX), disp%%spX, disp)
  disp<-ifelse(disp>spX, disp%%spX-spX, disp)
  popmatrix[,"X"]<-disp
  popmatrix
}

# Calculates Environmental Variance
VE<-function(h, VP)
  {
  (1-h)*VP
  }
##### #


sum.metrics<-function(popmatrix, spX, bw=10)
  {
  bins<-floor(min(popmatrix[,"X"])):ceiling(max(popmatrix[,"X"]))
  out<-.Call("sum_metrics", R_X=popmatrix[,"X"], R_H=popmatrix[,"H"],
             R_D=popmatrix[,"D"], R_n=nrow(popmatrix), R_bins=bins, R_nbins=length(bins),
    		 R_spX = spX, R_bw=bw)
  colnames(out)<-c("Nx", "MeanD", "MeanH", "SDD", "SDH")
  out<-cbind(X=bins, out)
  out
  }
##### #


ca.metrics<-function(popmatrix, spX)
  {
  bins<-floor(min(popmatrix[,"X"])):ceiling(max(popmatrix[,"X"]))
  cn<-c("Nx", "MeanD", "MeanH", "SDD", "SDH")
  out.sp<-.Call("sum_metrics", R_X=popmatrix[,"X"], R_H=popmatrix[,"HI"],
                R_D=popmatrix[,"D"], R_n=nrow(popmatrix), R_bins=bins, R_nbins=length(bins),
                R_spX=spX, R_bw=1)
  out.ind<-c()
  for (ii in 1:ncol(out.sp))
    {
    out.ind<-cbind(out.ind, approx(x=bins, y=out.sp[,ii], xout=popmatrix[,"X"])$y)
    }
  colnames(out.ind)<-cn
  matrix(out.ind, ncol=length(cn), dimnames=list(NULL, cn))
}
##### #

### Hassel-Commins population growth where b=1 (contest competition) #

hass.comm<-function(lambda, K, N)
  {
  lambda/(1+(lambda-1)*N/K)
  }
##### #






